subroutine shock_wave(u1,p1,rho1,u2,p2,rho2,D)
implicit none
!comput the velocity and rho

    real(8) D,p1,p2,u1,u2,rho1,rho2
    
    D=u1+(p2-p1)/(u2-u1)/rho1
    rho2=rho1*(D-u1)/(D-u2)

endsubroutine